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Abstract. In this paper we describe a new method for studying the hydrodynamical problem of a planet embedded in a 
gaseous disk. We use a finite volume method with an approximate Riemann solver (the Roe solver), together with a special way 
to integrate the source terms. This new source term integration scheme sheds new light on the Coriolis instability, and we show 
that our method does not suffer from this instability. The first results on flow structure and gap formation are presented, as well 
as accretion and migration rates. For Mp< 0.1 Mj and Mp> 1.0 Mj (Mj = Jupiter's mass) the accretion rates do not depend 
sensitively on numerical parameters, and we find that within the disk's lifetime a planet can grow to 3 — 4 Mj. In between these 
two limits numerics play a major role, leading to differences of more than 50% for different numerical parameters. Migration 
rates are not affected by numerics at all as long as the mass inside the Roche lobe is not considered. We can reproduce the 
Type I and Type II migration for low-mass and high-mass planets, respectively, and the fastest moving planet of 0.1 Mj has a 
migration time of only 2.0 10* yr. 
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1. Introduction 

, It is now generally believed that planets are formed out of 

■ the same nebula as their parent star. When this cloud of gas 
' collapses into a protostar, conservation of angular momen- 
tum leads to the formation of an accretion disk around the 

I star. These disks are inde ed observed around T-Tauri stars 

■ jBeckwith & Sargen? '1996^, and it is within these disks that 
I planet formation should take place. 

In standard theory, terrestrial planets as well as the rocky 
' cores of gas giant planets arise slowly from collisions of dust 
particles. When a protoplanet reaches a certain critical mass 
(w 15 M0), it can no longer sustain a hydrostatic atmosphere, 
and dynamical accretion sets in. Eventually, this accretion pro- 
cess forms a gaseous envelope around the core, of a mass com- 
parable to Jupiter. This scenario is known as the 'core accretion 
model'. 

Recen tlv. IBossI ( fl997| ) revisited a model originally pro- 
posed by ICameroij ( Il978h . the 'core collapse scenario', in 
which a giant planet is formed in much the same way as a star 
through a gravitational instability in the disk. However, Rafiko v 
( I2OO5I) pointed out several problems with this scenario, and it is 
yet undecided which of these two scenarios is coiTect. In both 
cases we end up with a newly formed giant planet still em- 
bedded in the protoplanetary disk. It is this interesting stage of 
planet formation that we focus on in this paper. 

Send ojfprint requests to: S. J. Paardekooper 
e-mail: paardeko@ strw . leidenuniv . nl 



The disk and the planet interact gravitationally with each 
other. The planet perturbs the disk through tidal forces, break- 
ing its axisymmetry and, if the planet is massive enough, opens 
up a gap in the disk (Lin & Papaloizou 1993). The perturbed 
disk in its turn exer ts a torque on the planet, lead ing possibly 
to orbital migration taol dreich & Tremainell983 ). Inward mi- 
gration is generally believed to be the mechanism responsible 
for creating 'Hot Jupiters' : giant planets, comparable to Jupiter, 
with very small semi-major axes (< 0.1 AU). 

The process of gap formation raises two interesting ques- 
tions. First: does the presence of the gap prevent further gas 
accretion? If it does, this puts serious constraints on the max- 
imum planet mass that can be reached. Secondly: how is the 
orbital evolution of the planet affected by the gap? If all plan- 
ets move inward on much shorter time scales than the typical 
lifetime of the disk, how come we still have Jupiter in our Solar 
system at 5.2 AU? 

The accr etion i s sue h as been investigated in deta il nu- 
merically by iKIevI ( 1 19991) and by iLubow et all (Il999l) . and 
they find that accretion can continue through the gap, allow- 
ing more massive planets to b e formed. Migration has been 
investigated both analy tically dOoIdreich & Tremainel Il98(i 
iLin & Papaloizou 1986), recently by 'Rafikov ('2002'), and nu- 
merically by Nelson et al. ( 2000) and Kiev (2000) . and it seems 
that migration is always inward, allowing for Hot Jupiters but 
not for 'Cold' J upiters. Recently Masset & Papaloizou 
and lArtvmowicz, (.2004) showed that there exists a runaway 
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regime in which the orbital radius of the planet evolves on very 
short timescales, with possible outward migration. 

However, the problem is complicated analytically, and hard 
to do numerica lly as well. A nice illustration of this is given 
bv lKlevl ( ll998l) . who showed that the use of a corotating coor- 
dinate frame can lead to non-physical evolution in some nu- 
merical codes. It therefore seems appropriate to introduce a 
numerical method which is new in this field of research, and 
which should handle rotating coordinate frames and the ef- 
fects of gravity in a better way. We use a finite-volume method 
based on an approximate Riemann solver for arbitrary coordi- 
nate frames, which can is specifically aimed to treat disconti- 
nuities in the flow correctly. Also, the numerical scheme is con- 
servative, meaning that the method conserves mass, momentum 
and angular momentum up to second order. 

In this paper, we aim to give a full description of RODEO 
(ROe solver for Disk-Embedded Objects), to show some dif- 
ferences in the gas flow in the disk and to present results on 
accretion and migration rates. The focus will be more on nu- 
merical effects than on physical effects, i.e. we do not consider 
different disk structures or diffe rent magnitudes of an a noma- 
lous v iscosity as was done by Brvden et al.l ( 1 19991) and iKlevI 
l ll999l) . Instead, we consider the effect of various numerical 
parameters to assess their importance regarding gap formation, 
accretion and migration rates. 

We focus on the two-dimensional, vertically integrated 
problem, which is formally only coiTect for planets for which 
the Roche lobe is larger than the disk scale height. For our 
disk, this comes down to Mp> Mj. However, because two- 
dimensional simulations are far less expensive in terms of com- 
puting time than three-dimensional simulations this allows us 
to do a detailed numerical study and compare our findings to 
previous two-dimensional results. 

We start in Sect. |2lby describing the physical model that 
was used. Next, we focus on the numerical method in some 
detail in Sect. |3] In Sect. |3 we show the results of some test 
problems and in Sect. |5] we give the initial and the boundary 
conditions. We give our results on gap formation, accretion and 
migration in Sects. |6l0and|8] Section|9]is reserved for a short 
summary and the conclusions. 

2. Physical Model 

2. 1. Basic equations 

Protoplanetary disks are fairly thin, i.e. the ratio of the vertical 
thickness H and the radial distance from the center r is smaller 
than unity. Typically H/r — 0.05. We can therefore vertically 
integrate the hydrodynamical equations, and work with verti- 
cally averaged state variables. We will use cylindrical coordi- 
nates (r, (/)), with the central star of 1 M0 located at r = 0. 

The flow of the gas is determined by the Euler equations, 
which express conservation of mass, momentum in all spatial 
directions, and energy. These equations can be written in a sim- 
ple form: 



Table 1. Definition of symbols 



dW dF dG 



r radial coordinate 

(j> azimuthal coordinate 

E surface density 

Vr radial velocity 

U0 angular velocity in the corotating frame 

p vertically integrated pressure 

"I> gravitational potential 

Q, angular velocity coordinate system 

Tp orbital radius of the planet 

P orbital period of the planet 



where W is called the state, F and G are the radial and the po- 
lar flux, respectively, and S is the source term. The state con- 
sists of the conserved quantities, and Eq. Q expresses the fact 
that any change in time of the state in a specific volume is due 
to flow across the boundary of the volume (the flux terms) or 
due to a source inside the volume (the source term). We will 
omit the energy equation for now, because we will work with 
a simple (isothermal) equation of state, which does not require 
an energy equation. The method still includes the energy equa- 
tion as an option, however, and this way we have another way 
of simulating isothermal flow: a run w ith a very low adiabati c 
exponent F, as was done for example in iNelson & Ben3 ( l2003l) . 
We will present the method only for the isothermal equations; 
for the full method in cluding the energy equation we refer to 
lEulderink & Mellemai d 19951) . 

In cylindrical coordinates, the vectors W , F, G and S can 
be written in the following form: 



W = r(E, Ewr, Sw^) 

F = r{Y.Vr, +p, SWrW^) 







(2) 
(3) 
(4) 



(5) 



The symbols used are listed in Table l2.ll 

Note that in this form there are no derivatives of depen- 
dent variables in the source term, keeping S finite even in the 
presence of discontinuities. Also, there are no viscous terms 
present, because we deal with viscosity in a separate way (see 
Sect.O- 

The gravitational potential contains terms due to the central 
star and due to the direct and indirect influence of the planet: 



Here r-i = 



r2 



(6) 



r-p I denotes the distance to the planet. Because 
this potential is singular at the position of the planet, and to 
account for the threedimensional structure of the disk, we use 
a smoothed version of r2 : 



dt 



dr 



(1) 



r2 



2r Tp cos( 



+ 



(7) 
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The smoothing parameter e is taken to be a certain fraction of 
the Roche lobe of the planet: 



e = b Rr =hrr, 



(8) 



Throughout this paper we have used h — 0.2. The results do 
not depend much on its value, as long as & < 1. The indirect 
term in the potential (the last term on the right-hand side in Eq. 
|6j arises due to the fact that a coordinate system centered on 
the central star is not an inertial frame, because the star feels 
the gravitational pull of the planet. A similar (small) term for 
the disk is omitted for simplicity. 

2.2. Equation of state 

Equation ([3 needs to be complemented by an equation of state. 
We assume an isothermal equation of state: 



(9) 



This choice is based on the assumption that the gas is able to 
radiate away all its excess energy very efficiently. In vertical 
hydrostatic equilibrium, the isothermal sound speed Cs is given 
by: 



H GMp 



© 



(10) 



2.3. Viscosity 



The nature of the viscosity in accretion disks was unknown 
for a long time. Currently, the Magneto-Rotational Instability 
(MRI) is the best candidate for providing a turbulent viscosity 
tBalbus & Hawlev 1990) . In regions where the ionisation frac- 
tion in the disk is high enough to sustain the MRI one would 
need to do Magneto-Hydrodynamics (MHD) to study planet- 
disk interaction, as was done in Nelson & Papaloizou (2003 ). 
MHD simulations are very expensive however, and it is not 
yet clear if the MRI operates throughout the disk. In so-called 
dead-zones cGamxnie 199 6) there may be no turbulent viscos- 
ity at all. However, in order to compare our results to previ- 
ous studies we include an anomalous turbulent viscosity in our 
models. 

Viscosity comes in as two extra source terms for the mo- 
mentum, one in the radial and one in the azimuthal direction. 
We deal with these source terms separately in the numerical 
method (see Sect.|3j. 

The f orm of these source terms can be found for example 
in lD'Angelo et al. (2002): 



/r 



1 d{rS„) ^ 1 d{Sr^) _ 
r dr r d(j) r 

r dr 86 



where the components of the viscous stress tensor S are: 



dr 



1 



(11) 
(12) 

(13) 





( dVrh 

^ [-di 


r 




/I r)ii 

\r d(f> 


+ r 

dr 


V • V 


1 d{rv^) 


dv^ 


r dr 


'86 



(14) 
(15) 
(16) 



W e can either use an g-presc ription for the viscosity parameter 
ly dShakura & SunvaevI 19731) : 

ly = aCsH (17) 

or we can take v constant throughout the disk. Note that for 
the a-disk with constant aspect ratio H/r, v oc -y/r, leading to 
enhanced viscosity in the outer disk. 

3. Numerical Method 

We can integrate Eq. Q over the finite volume of a grid cell to 
obtain: 



^— — if dr [ d6 (Ty"+^ - W") + 
\rA(j) J J 

J dtj d6 {F,+ i/2,j - -F^-l/2j) + 

j dt j dr (Gij+i/2 - Gjj-i/2) - 
dt dr d6S) ^0 



(18) 



Here Afj means the physical quantity A, evaluated at time 
index n, at coordinates The volume term (r^sin6' for 

spherical coordinates, r for cylindrical coordinates) of the grid 
cell is already present in W, F and S. A second order accurate 
integration scheme for this equation is: 



l.i{wr+^ - {WD 



1 



i+l/2,j ^ i-l/2,j) 
A(j) iJ + l/2 '^i,j-l/2) 



pn+l/2 



(19) 



where () denotes arithmetic mean. A numerical scheme like 
this is called conservative because the conserved quantities are 
indeed conserved by the numerical method: what goes into one 
cell must come out of another. 

We will use Eq. ( I19l l to find an update for the state. What 
is left to be done is to find expressions for the interface fluxes 
(Sect. ^nd to account for the source terms (Sect. ll!2t . 



3. 1. Roe solver 

First, we use the technique of dimensional splitting to obtain 
the two one-dimensional equations 

dW ^d£_^ 
8t dr 
8W 8G 

77- = Y (20) 



dt 



4 



Sijme-Jan Paardekooper and Garrelt Mellema: RODEO: a new method for planet-disk interaction 



The order in which these equations are solved is alternated to 
avoid systematic numerical effects (Strans 1968). We have the 
freedom of choosing the splitting of the source term any suit- 
able way. We discuss a special way in Sect. 13.21 

From now on, we focus on the radial direction; the az- 
imuthal integration is done in a similar way. We split Eq. (|20j 
once more to obtain an equation without source terms, and 
solve t his equation with a method originally prop osed by Roe 
(Il98lh . and extended by ' Eulderink & Mellemal (119951) to a 
general relativistic method. The non-relativistic limit of this 

method yields a solver for the Euler equations in arbitrary co- Fn — Fl = 6^6 
ordinate s. 



The interface flux is related to the interface state by; 
F = AW = AQC = QVC = ^ bkSk (26) 

k 

where hk — ^ka,k- We present the exact expressions for 
the eigenvalues, eigenvectors and projection coefficients in 
AppendixlAl 

We can project the flux difference at the interface we are 
considering on the eigenvectors of A: 

(27) 



3.1 .1 . Characteristics 

Given the state (or the flux) immediately left and right of an in- 
terface of two grid cells, the Roe solver computes the resulting 
interface flux by solving: 



dW dW 



dt 



dr 



(21) 



where ^(H^) — dF/dW is the Jacobian matrix. Roe's central 
idea is to approximate ^( W^) by a constant matrix A. Since the 
Euler equations are hyperbolic, this matrix A has eigenvectors 
efc and corresponding eigenvalues Afc. These can be used to 
diagonalize A: 



Q-^AQ = V 



(22) 



where Q is the matrix with right eigenvectors of A. Now we 
can cast Eq. i2l\ into characteristic form: 



9C dC „ 

at or 



where C is the vector of characteristic variables, defined by 
dC = Q-^ dW (24) 



and V is the diagonal matrix with eigenvalues of A. From ( I23> 
it is easy to see that dC = along the path r ~ Vt. These 
paths are called characteristics, and are essential in the study 
of hyperbolic differential equations. Discontinuities (shocks) 
travel along characteristics, and the domains of influence and 
dependence are bounded by the characteristics. 

We can integrate Eq. i23\ to find a relation between the 
state and the characteristic variables: 



W = QC = J2akek 



(25) 



where is the k-th characteristic variable. As the are used 
to project the state on the eigenvectors of A, they are called 
projection coefficients. 

We can use the projection coefficients for calculating the 
interface flux, because we know that they are constant along 
characteristics. If we can figure out from which point in space 
the characteristics that cross the interface originate at the cur- 
rent time, we can just calculate the projection coefficients at 
that point in space and use Eq. (I25> to find the state at the inter- 
face. 



where Fji and Fl aie the fluxes immediately right and left of 
the interface, respectively. Then the first order interface flux is 
approximated by: 

K+i/2 ^liFL+Fi,)-^J2 '^^^fc^fc (28) 
where ak — sign(Ak). 



3.1.2. Flux limiter 

The first order expression for the flux, Eq. J28> . assumes a jump 
in the projection coefficients. That is: depending on the sign of 
the eigenvalue Afc, we use the bk coiTesponding to either the left 
or the right state. This procedure is correct if there is a shock 
right at that interface, so the state makes a jump there. In re- 
gions of smooth flow however, linear interpolation between the 
different bk seems to be a better approach. To switch between 
the two kinds of interpolation a flux limiter tp is used. We fol- 
low the suggestion by Roe (1985), in which one uses the ratio 
of the state difference at the interface ao — [ak]i+i/2 to the 



(23) upwind state difference: 



(a/c)i_i/2 if Afc > 0; 
(afc)j+3/2 if Afc < 



Let 



i/'(ao, flu) — max(0, min(pau, max(ao, min(au,pao)))) + 

min(0, max(pa„, min(ao, max{au , pao)))) (29) 

This form of the flux limiter allows us to tune the way the 
method deals with sharp discontinuities through the parame- 
ter p. When p — 2, for example, the flux limiter is called 
'Superbee'. This limiter tends to create very sharp shocks, but 
tends to steepen shallow gradients, leading to numerical prob- 
lems such as under- and overshoots. For p = 1 the limiter is 
called 'minmod', which is more diffusive. To get the best of 
both worlds (sharp shocks but no under- or overshoots), we set 
p ~ 1.5 in all simulations. 

With this flux limiter, the second order interface flux be- 



^1/2 



Fr)-^ y^(o-fcQfc-(crfc-t^fc)i/'fc)Afcefc(30) 



where i/k — ^k^t/^r, the coefficient needed for linear inter- 
polation. When ipk = (shock), Eq. ( I28> is retrieved, while for 
i/'fe = 1 (smooth flow) the ak in Eq. (I28> is replaced by I'k and 
we interpolate linearly between the projection coefficients. 
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3.2. Stationary Extrapolation 

The usual way for dealing with source terms is to split Eq. |20| 
once more to end up with an ordinary differential equation for 
the state: 



dW 

~dt 



dF 

dr 



= 



thereby assuming that the flux is constant in space. The ordi- 
nary differential equation can be solved with any second-order 
integration scheme to yield a second-order update for the state. 
However, it is also possible to take the opposite approach: 



dW 
~dt 



= 0, 



dF 

dr 



= X 



(32) 



This is a method we call stationary extrapolation, because it 
assumes a stationary state within one grid cell. This equation is 
more difficult to solve, but it has certain advantages: 

- The Roe solver solves a Riemann problem at the interface, 
a configuration with two stationary states separated by a 
discontinuity. Using actual stationary states on either side 
of the interface ensures that the Roe solver deals with a 
genuine Riemann problem. 

- Physical stationary states are recognized. When the actual 
states are stationary, the Roe solver will produce no (un- 
wanted) numerical evolution of these states. This property 
is related to the numerical instability found by KlevI (Il998l) 
concerning the Coriolis force. We will see below that a 
scheme using stationary extrapolation does not suffer from 
this instability. 

If we split the source terms appropriately: 



— 2I]?;r(r2 + Vff,) 



Y = 




(34) 



the i ntegratio ns can be done analytically 
( lEulderin k & Mel lemal Il995l) . For isothermal, stationary 
flow in the radial direction Eq|32]can be rewritten: 



d_ 

dr 



■^vj) + $ + ( 
r^{v^ + n) 



■log(S) 



(35) 



That is: the mass flux, the Bernouilli constant and the spe- 
cific angular momentum are invariant. From these invariants 
the state at a cell interface can be computed from the state at the 
cell center. However, this procedure is computationally expen- 
sive, therefore it is feasible to adopt a first order approximation 
to calculate the interface fluxes: 



Fr-i/2M = F^- ^ArX, 



-1/2, i 



(36) 



This procedure can conserve station ary sta tes up to second or- 
der (Mellema et al. 1991). However. IklevI fl 998 ) showed fliat 
the angular momentum in disk simulations deserves special at- 
tention. We illustrate this with a simple example. Consider the 
radial transport of angular momentum: 



(31) d_ 

dt 



d 



(37) 



For stationary radial flow the mass flux and the specific angular 
momentum are constant: rSu^ — D and r^(w0 + 0) = L. 
For simplicity we assume here that D is constant for the whole 
computational domain. Consider the interface between cells i 
and i — 1. Stationary extrapolation from the cell centers to the 
interface i — 1/2 leads to the fluxes 



-1/2 



Ft, 



D{- 



j-1/2 



(38) 



(39) 



With Vj. > 0, the first order interface flux produced by the Roe 
solver is ^i_i/2 = Fl - Therefore 

„2 



-1/2 



rtu2{FLlD + ^l) = U^, 



(40) 



Stationary extrapolation from the interface back to the cell cen- 
ter i gives the final flux coming from the left for cell i: 



From the same analysis for interface i 
going to the right for cell i: 



(33) Fr^, ^D{- 



(41) 

1/2 we find the flux 
(42) 



The first order update for the state is: 
At 

Ar r; 



A(rS«^) = — (Fi,, 
Ar 



R,i) 



AtD 

^[L,^l 



(43) 



The change in angular momentum J 
this update is given by: 



AJ 



M D 



r^i;(w0 + Vl) due to 



(44) 



It is e asv to see that the conservative formulation of iKlevI 
(Il998h : 



(45) 



leads to the same change in angular momentum. This shows 
that this new method does not suffer from the numerical insta- 
biUty du e to the expli cit treatment of the CorioUs force that was 
noted by Kiev '(' 19981) . Note however that this does not hold for 
the approximate extrapolation of Eq. |36l Therefore we always 
use the exact extrapolation for the specific angular momentum, 
while keeping the method of approximate extrapolation to deal 
with the other source terms. This way we do not need to do the 
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computationally intensive full stationary extrapolation, while 
keeping angular momentum nicely conserved. 

Stationary extrapolation therefore provides a different point 
of view on the Coriolis instability: the failure of numerical hy- 
drodynamic schemes to properly conserve angular momentum 
can be traced back to the failure to recognize stationary s tates. 
Rewrit ing the angular momentum equation as done by iKlevI 
( Il998l) is a way to solve this problem, but it only deals with 
the Coriolis source term while similar problems may exist due 
to the other source terms as well. Therefore we feel that it is a 
good idea to integrate all source terms using stationary extrap- 
olation. 

One disadvantage of stationary extrapolation is that it is not 
known a-priori if the assumption of stationary flow is valid. 
This is especially important when the source terms are large, 
in our case very close to the planet. When dealing with more 
massive planets (> 0.5 Mj) the assumption of stationary flow 
along the coordinate axes breaks down, leading to numerical 
instabilities. A physical explanation is that the gas in this case 
will try to orbit the planet, perhaps forming a Keplerian disk, 
rather than to orbit the star. The flow is still semi-stationary, but 
only in a cylindrical coordinate frame centered on the planet, 
and this flow is at some points even orthogonal to the stationary 
flow in the frame of the star. Therefore, it seems appropriate to 
treat these larger source terms as external (see section l3.2.1> . 
This transition from stationary to non-stationary extrapolation 
is applied smoothly at a typical distance i?R from the planet 
with a sin^ ramp. 

3.2.1. External source terms 

The stationary extrapolation deals with the geometrical source 
terms (including gravity, which is merely a geometrical effect 
in general relativity). Any other source term which we will call 
Z is better integrated using Eq. i31\ : 



W"+^ = W" + AtZ 



(46) 



In our case, Z consists of the viscous source terms. The deriva- 
tives in these terms are calculated using a finite-difference 
method, and the resulting source is fed into Eq. (I46> . 

3.3. Adaptive Mesh Refinement 

Resolution is always an issue in numerical hydrodynamics, and 
a compromise has to be found between resolution and compu- 
tational cost. Adaptive Mesh Refinement (AMR) is a very eco- 
nomic way of refining your grid where the highest resolution is 
needed, whereas keeping large parts of the grid at low resolu- 
tion. 

We have i mplemented the PARAMESH algorithm 
(iMacNeice et alJEoQO) to be able to resolve the region near 
the planet accurately. Usually the refinement criterium is taken 
to be the second derivative of the density, but we can also 
predefine the region that has to be refined. When running in 
this mode, and switching off additional refining and derefin- 
ing, the algorithm works like the nested grid technique of 
[p'Angelo et al.. (.2002) . However, since our grid can be fully 




Fig. 1. The evolution of the axisymmetric gas ring, initially lo- 
cated at x=l, shown at r = 0.05, and after 50 and 100 or- 
bits. For the last two the numerical results are also plotted (di- 
amonds and triangles, respectively. 



adaptive it is suited to let the planet migrate through the disk 
while keeping a high resulution near the planet. This will be 
the subject of a future paper. 

We use linear interpolation to communicate boundary cells 
between the different levels of refinement. The implementation 
of the monotori i sed ha rmonic mean ( van Leer 1977) used by 
iD'Angelo etaT]ll2002h made no significant difference. For ev- 
ery boundary between different levels we reset the flux into the 
lower level in order to conserve mass and momentum across 
the boundary. 

3.4. Accretion 

We foflow b' Angelo et alJ ll2002l) in modeling the accretion by 
taking away some fraction of the density within a distance of 
Tacc of the planet. Basically, the density is reduced with a factor 



1 - / At, 



(47) 



where At is the magnitude of the time step. The two parameters 
that determine the accretion rate are race and /, where / is 
taken to be twice its value in the inner half of race- The mass 
taken out is monitored, assuming it has been accreted onto the 
planet, but without changing the dynamical mass of the planet. 

Equation ^\ can be seen as a first order approximation to 
the solution of the differential equation 



dp 



p 



with the solution 



p{t + At) - p(t)e-^*/^- « p{t) (1 - At/Tacc) 



(48) 



(49) 



From this equation we see that the typical time scale race for 
emptying the accretion region is /^^. 
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4. Test Problems 

We tested the RODEO method against several hydrodynami- 
cal t est problems. First of all, the standard Sod shocktube (Sod 
11978 ). in one and in two dimensions. For the two-dimensional 
problem, we have placed the initial discontinuity diagonal 
along the grid, making it a genuine two-dimensional prob- 
lem. The analytic solution was nicely recovered, except at the 
boundaries where some reflections were seen. 

A more complicated te st is the windtunnel with step 
llWoodward & Colella 1984*). which was previously used as a 
testproblem by Mellema et al] (Il99ll) using the same method. 
The results agreed very well. Note that both tests do use an 
energy equation. 

4.1. Viscous ring spreading 

To test the implementation of the viscous terms, we set up an 
axisymmetric ring of gas in Keplerian rotation about a cen- 
tral object. If we neglect pressure forces, the evolution of the 
surface density, which initially is a delta function, is given by 
[Pringle(1981): 

, , 1 1+.^ f2x\ 

where x denotes the normalized radial coordinate r/ro, with 
the ring initially located at r = rg, and r denotes the dimen- 
sionless time (t = t/ty) in units of the viscous spreading time 
ty = rQ/12v. Ii/4 denotes the Bessel function of the first kind 
of fractional order 1/4. 

As the analytical solution neglects any forces due to a pres- 
sure gradient, we have to make sure that we set the temperature 
to a very low value, or else pressure waves will dominate the 
pattern. We have used the same resolution as for simulations 
with the planet, and with a viscosity parameter a = 0.004 at 
a; = 1, comparable to the value used in the disk-planet simu- 
lations. We have placed the ring in the middle of the compu- 
tational domain of x G [0.27, 1, 73], with t = 0.05 initially 
because we cannot model a delta function numerically. 

We have also performed a test run of this problem using 
a = to check for numerical viscosity, and it turned out 
that this (unwanted) numerical diffusion was very much lower 
than any physical viscosity, at a low resolution of 128x384. 
Note that this also shows that we can maintain a stable disk in 
Keplerian rotation. This is a nice result, showing that our con- 
servative scheme conserves angular momentum very well. 

Figure [0 shows the results for a = 0.004. It is clear that 
the numerical solution agrees very well with Eq. (I50> . indicat- 
ing that the implementation of the viscosity works. Only at the 
boundaries there are minor deviations from the analytical solu- 
tion. 

5. Model design 

Throughout we use non-dimensional units. The unit of distance 
is the orbital radius of the planet, and the unit of time is the 
inverse orbital frequency. Note that this differs by a factor 2tt 
from the orbital period of the planet. 



5.1. Initial Conditions 

The base resolution in all our simulations is {n,:,n^) = 
(256, 768), which leads to square cells near the position of the 
planet. We set h = H/r = 0.05, which determines the temper- 
ature through Eq.[roi We take a constant initial surface density. 
Because we do not consider the self-gravity of the disk the den- 
sity can be given in any desirable unit, and we normalize it to 
1 at the planet's position. We set the radial velocity to zero, al- 
though another possibility is to take an exact a-disk initially. 
The initial angular velocity is set to the value for a Keplerian 
disk, with a small connection for the pressure force. 

We take a constant kinematic viscosity i/ = 10^^, corre- 
sponding to a = 0.004 at the location of the planet. We put in 
a planet at location (r, </)) = (1, ti")- This planet can be allowed 
to accrete matter (without changing its dynamical mass). 

5.2. Boundary conditions 

The inner and the outer boundary are located at 0.4 and 2.5, re- 
spectively. Imposing proper boundary conditions is not trivial, 
because waves are continuously hitting both edges of the grid, 
changing the sign of the radial velocity. A standard outflow 
boundary assumes that the velocity is always directed outward, 
and in combination with outgoing waves this leads to insta- 
bilities near the inner radius of the disk. A reflecting boundary, 
on the other hand, will lead to reflected waves traveling into the 
computational domain, which interact with the outgoing waves, 
destroying the flow structure. 

Therefore we decided to take a non-reflecting boundary, 
following Godon ( 1996). This boundary treatment is based on 
characteristics, and as a result the implementation is particu- 
larly simple in our numerical scheme. The idea is to impose 
the boundary conditions on the characteristic v ariables , rather 
than on density, momentum or energy directly. iGodorJ (Il996l) 
showed that this leads to a non-reflective boundary. 

A ghost cell has to be updated according to incoming char- 
acteristics: when a characteristic leaves the last regular cell into 
the ghost cell, the corresponding characteristic variable of the 
ghost cell is updated accordingly. When a characteristic enters 
the ghost cell from outside the numerical domain, the charac- 
teristic variable is updated using the unperturbed (Keplerian) 
value. Using the Roe solver, this comes down to having two 
ghost cells: the outer one is never updated, so density and ve- 
locities remain at the initial values, and the inner one serves as 
the actual ghost cell. The latter is automatically updated by the 
Roe solver us ing characteris tic variables, thereby following the 
suggestion by 'Godon' 7l996V 

It turns out that these boundary conditions work especially 
well with stationary extrapolation, because then the hydrody- 
namics part as well as the boundary prescription are deal- 
ing with small fluctuations over a stationary background state. 
When we switch to the ordinary integration of source terms of 
Eq. lathis is no longer true. In this case small numerical arte- 
facts can be seen near the boundaries. 

Therefor e, as an addition we ca n add a wave-damping 
mechanism jPe Val-Borro et al1l2005l) that operates in the re- 
gions 0.4 < r < 0.5 and 2.1 < r < 2.5. In these regions the 




state is relaxed to the initial (stationary) state on a time scale 
varying from infinity at r = 0.5 and r = 2.1 to 1 orbital period 
at the location of the boundary. Using this prescription, outgo- 
ing waves are damped gently as they approach the boundary of 
the computational domain. 

6. Gap Formation 

Only massive planets are able to open gaps in gas disks. There 
are two criteria that have to be fullfilled: firstly, the torques 
arising due to the presence of the planet must be able to 
overcome the visco us torques, leading to a minimum mass of 
(lBrvdenetal.ll999l) : 



Mr. 



(51) 



Secondlv. lLin & Papaloizoul(ll993b suggested that at a distance 
of Rr from the planet the planet's gravity should be more im- 
portant than pressure, which leads to another minimum mass 
for gap opening: 



3 h-^ M, 



o 



(52) 



For our disk, with v ~ 10^^ and h — 0.05 both criteria 
yield approximately the same minimum mass, namely 0.4 Mj. 
These criteria were shown to be valid within a factor of 2 in 
Ifirvden et al. ( 1999). 

In our simulations, we found the density reduction factor 
near the orbit of the planet to be more than 100 for a 1 Mj 
planet, 10 for a 0.5 Mj planet and 2 for a 0.1 Mj planet. All 



factors were measured after 500 orbits of the planet. In view 
of these results, we can confirm that both criteria give a good 
estimate of the minimum mas s for gap openin g. This seems to 
contradict recent results from lRafikovl ll2002h . who found that 
it takes only a fraction of Mmin.h to open a gap. This fraction 
depends on viscosity and disk mass, and it can be as low as 
0.1. However, in that analysis feedback from migration of the 
planet plays a role, while our planet remains fixed in the grid. 

Figure 12] (left panel) shows a greyscale plot of the surface 
density after 10 orbits of a 1 Mj planet, and all the basic fea- 
tures in the flow are already visible. We can see two trailing 
spiral arms leaving the planet: the inner arm moves all the way 
down to the inner boundary, while the outer arm reaches to 
about r = 1.5. Note also the secondary waves excited in the 
inner and in the outer disk. 

All spiral waves are stationary in the corotating frame. For 
this run we did not take out any matter inside the Roche lobe, 
so the density reaches very high values near the planet. 

The formation of the gap is quite a violent proces. There 
are structures visible in the corotating region and in the outer 
edge of the gap. The physical viscosity is able to damp these 
fluctuations, leading to a stable gap after 200 orbits (Fig. |2j 
right panel). 

The left panel of Fig. |3] shows the azimuthally averaged 
surface density for the same run at different times. The gap 
that is formed is approximately 0.5 rp wide, with two density 
bumps on either side. After 100 orbits the inner disk is starting 
to be accreted onto the central star. At all times the position of 
the planet is visible as a small spike at r = 1. 
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Fig. 3. Gap formation for a 1 Mj planet. Left panel: azimuthally averaged surface density for the case of exact stationary extrap- 
olation. Right panel: azimuthally averaged surface density for the case of approximate extrapolation. The indicated times are in 
planetary orbits. 
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Fig. 4. Azimuthally averaged density profile after 20 orbits for 
different implementations of the source term integration: ex- 
act stationary extrapolation, approximate extrapolation and no 
stationary extrapolation. 



The proces of gap formation i s numerically challenging, as 
was already shown bv lKlevllll998l) . In the remainder of this sec- 
tion we study the formation of gaps as a function of numerical 
parameters. 



6. 1. Source term integration 

First of all we focus on the integration of source terms. We con- 
sider four different methods. First of all our standard scheme, 
in which we solve for the angular momentum exactly, and inte- 
grate all other source terms using Eq.|36l Secondly, we have the 
approximate scheme, where all source terms are integrated us- 
ing Eq.m Next, we consider the case of no stationary extrapo- 
lation at all (see Ea.l31>. which is what ordinary Riemann-type 
schemes would use. Finally, we again integrate the angular mo- 
mentum exactly, and all other source terms using Eq.|2] This 
mimics the fix found bv lKlevI lll998h . and we will refer to this 
scheme as hybrid, because it combines the two extremes of sta- 
tionary extrapolation and no stationary extrapolation. 

In Fig. |3l we show the formation of the gap for the exact 
stationary extrapolation (left panel) and the approximate sta- 
tionary extrapolation (right panel). The difference between the 
two plots is quite dramatic. While both methods keep angular 
momentum conserved up to second order, approximate extrap- 
olation leads to a much wider and deeper gap. After 200 orbits 
of the planet, even the region 1.5 < r < 2.1 is participating in 
gap formation. From Fig.|3lit is clear that approximate extrap- 
olation, while it preserves stationary states up to second order, 
is not the way to go in case of a strongly rotating fluid. Only 
the exact stationar y extrapolation produces re sults comparable 
to other codes (seel Pe Val-Borro et aP ( 120051) '). 

For these runs the wave damping boundary regions were 
employed to be able to compare to runs with no stationay ex- 
trapolation at all. The effect of the damping zones is clearly 
visible after 200 orbits, especially at the inner boundary. 

Both models of Fig.|3]were also run in an inertial frame for 
comparison. For the case of exact extrapolation it made no dif- 
ference at all, while only small changes were seen for approxi- 
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Fig. 5. Azimuthally averaged density profile for a 1 Mj planet 
after 100 orbits for different flux limiters: p — 1.5 (standard 
value), p = 2.0 and p = 1.0 (see Eq.l29> 
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Fig. 6. Azimuthally averaged density profile for a 0. 1 Mj planet 
after 200 orbits in a non-viscous disk for different flux limiters: 
p = 1.5 (standard value), p = 2.0 and p = 1.0 (see Ea.l29l 

mate extrapolation. However, even in an inertial frame there is 
a Coriolis force term, see Eq.|5] so an inertial frame calculation 
is in th is case not a valid check for the instability found by Kleyl 
lll998l) . 

In Fig.0]we compare the azimuthally averaged surface den- 
sity profile for the four methods of source integration: standard, 
approximate, no stationary extrapolation at all, and the hybrid 
form. This figure shows that both the approximate extrapola- 
tion and the case of no stationary extrapolation overestimate 



the angular momentum transport, which leads to spurious gap 
widening. It is remarkable, however, that approximate extrap- 
olation does the worst job, worse even than the run without 
stationary extrapolation. This is especially clear for the inner 
disk (r < 1). It is interesting to no t e that in this case trying 
to deal with the issue raised bv lKlevI ( Il998l) in an approximate 
way actually makes matters worse than for the standard source 
term integration. 

The hybrid extrapolation reproduces the correct gap width 
and gap depth after 20 orbits, and follows closely the exact 
extrapolation almost everywhere. Only in regions where the 
strongest waves exist, near the planet and in the inner disk, the 
two methods give different results. Qualitatively hybrid extrap- 
olation gives similar results as a more diffusive flux limiter (see 
Sec.lO. 

We also performed the same simulations at a lower resolu- 
tion. For the case of exact and hybrid extrapolation no differ- 
ences were found, but for both other cases the width and struc- 
ture of the gap changed significantly, especially for the approxi- 
mate stationary extrapolation. It is clear from Eq.|36lthat the er- 
ror made in this approximate scheme will depend on the spatial 
resolution of the grid. On a grid of only {n^, n^) = (128, 384) 
the effect of gap widening is even more dramatic than seen in 
Fig.E] The inner disk is cleared away even faster, leading to an 
unstable situation near the boundary. For the case of no station- 
ary extrapolation the effect is not that severe, but still the gap is 
wider than in Fig.|4] 

6.2. Flux limiter 

The flux limiter is basically a switch between using a first or 
second order interface flux for the state update. Near shocks it 
should be first order, in smooth flow it should be second or- 
der. Applying a second order flux near a discontinuity leads to 
numerical smearing of the state, and therefore the flux limiter 
that switches the fastest to first order fluxes gives the least nu- 
merical diffusion. Here we study the effect of this numerical 
diffusion on the formation and appearance of the gap. 

In Fig.|5lwe show the gap structure for three different lim- 
iters (see Eq.l29li: p ~ 1.0 (minmod, most diffusive), p — 2.0 
(superbee, least diffusive) and p — l.h which we call soft su- 
perbee. In the outer disk all limiters give more or less the same 
result. But in regions where the waves induced by the planet 
are the strongest, the inner disk and close to the planet, we 
see clear differences. The diffusive minmod limiter damps the 
ingoing waves more than the other limiters, leading to an en- 
hanced surface density near r — 0.7. Close to the planet this 
higher diffusion also leads to an enhanced surface density, seen 
as a spike at r = 1 in Fig.|5l 

There are no large differences in the results obtained with 
the superbee limiter and its softer version. This indicates that 
our choice of p = 1.5 (see Sec. l3.1.2t is a good trade-off be- 
tween the two extremes minmod and superbee. It inherits the 
low numerical diffusion from superbee, as is clear from Fig.|5j 
while it is more stable against numerical overshoots. 

Since the different limiters imply different numerical vis- 
cosities, it is interesting to study how they influence the result 
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Fig. 7. Accretion onto a 1 Mj planet. Left panel: total accreted mass in units of the initial disk mass. Right panel: accretion rate, 
in units of disk masses per orbit. 



in the case when there is no physical viscosity added to the 
simulations. Then the situation changes drastically, because nu- 
merical viscosity starts to play a major role, in particular the 
length over which shocks are damped. As an illustration, we 
show the azimuthally averaged surface density for a 0.1 Mj 
planet in a viscosity-free disk in Fig.|6l for the three different 
flux limiters. It is clear that for the diffuse minmod limiter the 
waves are damped much faster, and the angular momentum is 
deposited much closer to the planet. For the superbee limiter, 
as well as its softer version, a much wider gap results. 

7. Accretion rates 

Now we turn to the problem of accretion onto the planet. The 
growth rate of the planet is important because it determines the 
ultimate mass the planet will reach. Two-dimensional studies of 
[p'Angelo et al. (2002) and Lubow et al. ( 1999) showed that a 
planet of 1 Mj grows approximately at a rate of 10~^ MdP~^, 
where Md is the disk mass within the computational domain. 
However, the study by Kiev ( 1999) done at lower resolution 
indicated an accretion rate more than ten times higher. 

In this section we look at three different mass regimes: 
high-mass planets, which open clear gaps in the disk, low-mass 
planets, which do not open gaps, and intermediate mass plan- 
ets, which create only small density dips around their orbit. 

7. 1. High mass planet 

We start by discussing the results for a 1 Mj planet. Because 
this planet is able to open up a gap in the disk, the accretion 
rate is determined by the amount of mass that flows through 
the gap ( Kiev 1999.) , and less by the density structure close to 
the planet. 

In Fig. we show the accreted mass and the accretion rate 
as a function of time for our standard resolution, accretion pa- 
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Fig. 8. Close-up on the logarithm of the surface density near 
the planet. Overplotted is the AMR mesh structure, where each 
block consists of 8x8 cells. 

rameters, source term integration and flux limiter (solid lines). 
Because we do not start with an initial gap the accretion rate 
is very high at the beginning of the simulation, dropping about 
one order of magnitude in 200 orbits to a value of 2 10^^ (disk 
masses per orbit). We ran one model to 500 orbits, and the fi- 
nal accretion rate turned out to be 1.5 lO^''. The fact that the 
accretion rate approaches a constant value after about 500 or- 
bits and the actual va lue found agree within a factor 2 with 
b'Angeloet 10 ( 120021) and lLubow etaP ( Il999l) . The planet ac- 
cretes approximately 10 percent of the total disk mass during 
the simulation. 
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Fig. 9. Accretion rates onto a 1 M j planet for three different 
methods of source term integration. 
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Fig. 10. Accretion rates onto a 1 Mj planet for three different 
flux limiters. 



7.1 .1 . Accretion parameters 

The accretion procedure is described by two parameters: the 
radius within which we take out material (race) and the value 
of / (see Eq.ETli. We can vary these to see if this influences the 
final accretion rate. 

In Fig. the accretion rates onto a 1 Mj planet are shown 
for the standard case (r'acc = 0.5 Rr and / = 0.5) and for 
parameters race = 0.1 Rr and / ~ 5/3 (t he curve lab eled with 
'higher f '). The standard set was used by Kiev ('199^, and the 
other case bv.D'Angelo et al. ( 2002 ). Note that the accretion ar- 
eas differ by a factor 25, while / varies only by a factor 3.3. So 
for identical density distributions close to the planet the stan- 
dard parameters yield 7.5 times more accreted mass during one 



time step. Despite this the final accretion rate is the same for 
both sets of parameters. 

Because of the different accretion radii any difference in 
accretion rate would imply that the flow within the Roche lobe 
is important for the accretion proces. Disk material that makes 
it to 0.5 Rr is accreted for the standard parameters, while it has 
to make it all the way to 0.1 Rr to be accreted in the second 
case. The fact that we do not see any differences indicates that 
the accretion rate is determined by the amount of mass the disk 
is able to supply, independent of local processes near the planet. 

7.1 .2. Resolution 

Our base resolution corresponds to approximately 8 cells per 
Rr, which means that the standard accretion area is resolved 
by only a few grid cells. This might be of influence on the in- 
ferred accretion rate, and therefore we performed simulations 
on different resolutions. 

First of all we lower ed the resolution with a factor of 2, 
the same resolution as in iKlevI lfl999.) . In this case we do not 
resolve the Roche lobe, so if the flow close to the planet is im- 
portant for accretion we would expect differences. However, 
since the two sets of accretion parameters did not yield differ- 
ent accretion rates we expect no big resolution effects. Indeed, 
Fig. shows that for this low resolution case (dotted line) the 
accretion rate is the same as for our standard resolution. This 
shows that resolving the flow within the Roche lobe is not im- 
portant for gap-opening planets, at least not for determining 
accretion rates. 

It is computationally very expensive to go up a factor 
of 2 in resolution on the whole grid, therefore we used our 
AMR module to refine the region around the planet. Figure |8] 
shows a close-up on the density pattern after 20 orbital periods. 
Overplotted are the Roche lobe and the grid structure. Each 
block represents 8x8 cells, so that we have approximately 1500 
cells within the Roche lobe. Therefore we can really resolve the 
flow inside the Roche lobe with this resolution. Nevertheless 
also this resolution yielded the same accretion rate of 2 lO^** 
MdP-i. 

7.1.3. Equation of state 

It is interesting to compare accretion rates for a truely isother- 
mal simulation and a run that does include an energy equa- 
tion, but at a very low adiabatic exponent F. This has been 
done before to mimic isoth ermal flow for planet-disk interac- 
tion ("Nels on & Ben3l2003l) . The basic idea is that for a low 
value of r the gas can be compressed without a large change 
in temperature. In order to check the validity of this approach 
regarding planet-disk interaction, we ran simulations including 
an energy equation but with a low value of T. 

We performed two simulations including an energy equa- 
tion, one with F = 1.00 1 and one with F = 1.01. The latter 
value was also used bv Nelson & Benj ( l2003l) . We found that 
the basic flow structures remain the same, and that the temper- 
ature profile does not change anywhere but very close to the 
planet. For F = 1.001 the temperature rises already by a factor 
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Fig. 11. Accretion onto a 1/64 M,j planet. Left panel: total accreted mass in units of the initial disk mass. Right panel: accretion 
rate, in units of disk masses per orbit. 



10, and for F = 1.01 by a factor of 60. This steep temperature 
gradient slows down the gas flow towards the planet consider- 
ably. 

Due to the higher temperatures close to the planet the ac- 
cretion slows down. Already for F = 1.001 the accretion rate 
drops by a factor of 2, and even a factor of 10 for F = 1.01. 
This shows that the accretion rate depends very much on the 
temperature near the planet, and that 'nearly' isothermal sim- 
ulations can produce very different results from truely isother- 
mal runs. 

This effect is different from the one described in iKlevI 
(Il999l) . who found that a polytropic equation of state leads to a 
reduction in the accretion rate. Because the temperature is pro- 
portional to the density in that case (for F = 2), the gap region 
is much cooler and therefore the viscosity is reduced when the 
a-formalism is used. In our case, there is a temperature rise 
very close to the planet, which leads to a pressure barrier that 
is able to slow down accretion considerably, even when F is as 
low as 1.001. Therefore we conclude that these kinds of sim- 
ulations are not able to mimic isothermal flow for this specific 
case, due to the deep potential well of the planet. 

7.1 .4. Source term integration 

In Fig. |3l we demonstrated that the way source terms are in- 
cluded has dramatic effects on the proces of gap formation. In 
this section we show that this also affects the accretion onto the 
planet. 

Figure|9lshows the accretion rates for three different meth- 
ods of source term integration: exact, hybrid and approximate 
extrapolation. See Sec.|6lfor their definition. Again, as in Fig. 
|3l approximate extrapolation is the odd one out, yielding a 4 
times lower accretion rate. This clearly has to do with the dif- 
ference in gap formation time scale, and again it is clear that 
approximate extrapolation is not the way to go. 



Both the exact and hybrid method give the same results. 
This is not surprising, because the region where the methods 
differ most is where the source terms are large: close to the 
planet. And as we mentioned before, this region is not impor- 
tant for the accretion rate. 

7.1 .5. Flux Limiter 

The flux limiter did not play a major role in gap formation (see 
Fig.O, only at the inner parts of the disk some differences can 
be seen. Figure^] shows the accretion rates for the three dif- 
ferent limiters discussed in Sec.|6l It is clear that the accretion 
rate does not depend at all on the flux limiter 

This can be understood if one realizes that it is the flow 
from the out er disk to th e inner disk that governs the accre- 
tion rate (see lKlevI lll999l) V But in the outer disk the waves are 
weaker than in the inner disk, so a different flux limiter should 
not change the mass flux from the outer disk very much. 

7.2. Low mass planet 

We now move to the other side of the planetary mass spectrum 
to investigate accretion onto planets that do not open gaps in 
the disk. Specifically we focus on a planet with mass Mj / 64. 

1/3 

Because Rr cx Mp the Roche lobe of this planet is 4 times 
as small as the Roche lobe of a jupiter-mass planet. 

In Fig.^Jthe solid line gives the accretion rate for our stan- 
dard parameters and 4 levels of AMR. For our base resolution 
the Roche lobe would only be resolved by only 1 grid cell, 
clearly not enough to study accretion. Figure ^3 dashed line, 
shows that 2 levels of refinement, yielding 4 cells per Rr, is 
still not enough to reproduce the result for 4 levels of AMR. 
The accretion rates for 3 and 4 levels of AMR agree very well, 
showing that at least 8 cells per Rr are needed for accurately 
modeling accretion. 
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Fig. 12. Accretion onto a 1/8 Mj planet. Left panel; total accreted mass in units of the initial disk mass. Right panel; accretion 
rate, in units of disk masses per orbit. 



All accretion rates have reached their final value after 30 
orbits. The model with 4 levels of refinement was run until 
200 orbits with no change in the accretion rate. This is because 
the planet does not open up a gap, which would take about 
200 orbits (see Sec.gJ. The value of 3.5 10"^ that we find i s 
in good agreement with the results of IP'Angelo et al.l ll2002l) . 
Note however that they need 7 grid levels, corresponding to 
approximately 6 levels of AMR for our simulations, while we 
need only 3 levels to obtain the same result. 

Because we have already discarded the approximate extrap- 
olation scheme in Sec. |6l and Sec. 16.11 we looked only at the 
difference between exact extrapolation and hybrid extrapola- 
tion for the low-mass case. However, we point out that in this 
case approximate extrapolation gave identical results. This is 
because the waves from this planet are too weak to make a dif- 
ference in angular momentum flow. From Fig. we see that 
both methods we considered for source term integration give 
identical results. In this case the planetary gravitational source 
terms are too small to cause a large difference in accretion rate. 

Also different flux limiters gave identical results. In Fig. 
II llwe only show the superbee-result, but the minmod limiter 
yielded exactly the same accretion rate. This was to be ex- 
pected, because in the limit of smooth flow (or, equivalently, 
no strong waves, as is the case for a low-mass planet) all lim- 
iters produce the same flux (see Eg. 1291 with ao = a„). 



7.3. Intermediate mass planet 

The masses that lie in between the two extremes we have con- 
sidered so far are perhaps the most interesting. In this regime 
the transition from linear disk response to gap formation takes 
place, and t herefore also the transi tion from Type 1 to Type 
II migration. Id ' An gelo et a l] |l2002l) find the highest accretion 
and migration rates here. Also, the cores of giant planets may 
well be in this mass range. 
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Fig. 13. Accretion rates onto a 1/8 Mj planet for the standard 
case (solid line), hybrid extrapolation (dotted line) and the min- 
mod limiter (dashed line). 



For this section we focus on a planet of mass 1/8 Mj, 
which makes its Roche lobe twice as small as for a Jupiter-sized 
planet. In Fig. [21 we compare the accretion rates for different 
resolutions. First of all, note that we need a relatively high res- 
olution to obtain convergence (3 levels of refinement, which 
amounts to 16 cells per Rr). This is already an indication that 
interesting things are going on. Secondly, we find an accretion 
rate that is about twice as low as was found bv iD' Angelo et all 
(2002). In view of the good agreement for the low-mass planet 
this is remarkable. 

In Fig. El we compare the standard model with the method 
of hybrid extrapolation and the minmod flux limiter. While for 
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Fig. 14. Left panel: accretion as a function of planet mass. Right panel: growth time (see Eq. I53t . Filled circles indicate our 
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iD'Angelo et al.. (.200Z) is shown by the dashed line. 



the high-mass planet as well as the low-mass planet these three 
methods gave identical results, for a 1/8 M j planet they differ 
significantly, giving a 66 % and 100 % higher accretion rate 
for the hybrid extrapolation and minmod limiter, respectively. 
With the most diffusive flux limiter we can reproduce the result 
ofE Angelo et aL(.2002.) . This shows that numerical diffusion 
is a very important issue in these kinds of simulations. 

It is not surprising that the diffusive minmod limiter in- 
creases the accretion rate, because it tends to smear out the 
strong density and velocity gradients near the planet, allowing 
mass to diffuse into the Roche lobe. For this planet local con- 
ditions determine the accretion rate, unlike the previous cases. 
Low-mass planets do not excite strong enough waves to alter 
their environment significantly, while high-mass planets open 
up gaps, and the accretion rate is therefore determined by the 
global evolution of the disk. Planets of intermediate mass, how- 
ever, do not clear a gap while they excite reasonably strong 
waves, making the dynamics very interesting. 

A similar story applies to the source term integration. We 
have seen no differences between exact and hybrid extrapola- 
tion for high and low-mass planets, while for our 1/8 M j planet 
the difference is quite dramatic. Again, the gravitational forces 
due to low-mass planets are too weak to cause a difference, 
while for high-mass planets the local conditions are irrelevant. 

7.4. Dependence on planetary mass 

In the previous sections we have looked in detail at three char- 
acteristic planetary masses. Here we focus on how the accretion 
rate depends on the mass of the planet. We have run simulations 
for planets from 0.5 up to 8 Mj, or, in other words, from 
deep in the linear regime to well above the gap-opening mass. 

In the left panel of Fig. [21 we show the accretion rate 
for all planetary masses, as well as the relation found by 



b' Angelo et all (l2002h . As was mentioned in Sec. 17.21 the re- 
sul ts for the low-mass plan ets agree very well with those found 

byE 

Angelo et al However, as soon as the disk re- 

sponse to the planet approaches the non-linear regime around 
Mp « 0.1 Mj the results start to differ significantly. At first, 
the accretion rate stays low, but around Mp w 0.2 Mj there is 
a strong rise leading to a maximum at 0.5 Mj, followed by a 
steep decline. 

The general features of the plot are consistent with the re- 
sults of ID' Angelo et al the accretion rate rises with 
planetary mass, with a maximum around 0.5 Mj, followed by 
a steeper decline. However, the rise as well as the decline are 
much more dramatic in our case, leading to a higher maximum 
accretion rate and lower accretion rates for the highest mass 
planets. 

The diamonds in Fig.ll4lrepresent the results obtained with 
the diffusive minmod limiter wherever they differed signifi- 
cantly from the standard case. We conclude that a diffusive flux 
limiter always tends to increase accretion onto the planet, espe- 
cially in the mass regime in which the transition from linear to 
non-linear disk response takes place. 

The difference between the results obtained with the min- 
mod flux limiter and our standard limiter can be as large as 
50% (see Fig. I14> . Because these differences are caused by 
purely numerical effects, one can interpret these as an esti- 
mate of the error in the values of the accretion rates in this 
mass regime. Keeping these eiTor estimates in mind, we see 
that our results are roughly in agreement with the relation 
found bv iD' Angelo et al. (2002). Note, however, that our stan- 
dard results always represent the lowest possible accretion rate. 
Different flux limiters or different source term integrations al- 
ways lead to a higher accretion rate. 
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calculations. 



for these calculations. 



The growth time scale is defined as the time it takes for a 
planet to double its mass at a given accretion rate: 

^ Afp 
- dMjdt 

In order to express Tg in years we took a disk mass of 0.0035 
Mq. When for a given mass Tg becomes comparable to the total 
disk life time the planet can not grow beyond this mass. This 
limiti ng mass is of the order of several Mj jP' Angelo et alJ 
I2OO2I) . However, they do not consider planets more massive 
than 1 Mj. 

In the right panel of Fig. we show Tg for all 
planet masse s. Again, the dashed curve shows the fit from 
Id 'Angelo et al. (2002 ). Because of the logarithmic scale on 
the y-axis the differences are much less pronounced. However, 
we still see the two different regimes and the sharp transition. 
Very important is the steep rise in Tg for the highest masses. 
Assuming that the planet spends maximal 10^ years inside this 
disk the maximum mass it can reach is 4 Mj. Note that if we 
exclude the planets more massive than 1 Mj the slope of the 
growth ti me scale is in excellen t agreement with the relation 
found by 'D'Angelo et al.' ("2002^. Extrapolating this slope we 
would find that the maximum mass that can be reached is about 
10 Mj. This shows that it is really necessary to include the 
highest mass planets to obtain a reliable estimate of the max- 
imum planetary mass that can be reached through disk accre- 
tion. 

The decline in accretion rate that we find for high-mass 
planets is steeper than found by Lubow et al, (1999) . who used 
the ZEUS code for their simulations. We have shown that this 
is not due to the flux limiter or the source term integration, 
and the reason for this difference may be found in the intrinsic 
difference between the finite-difference approach and Riemann 
solvers. However, this is impossible to verify within our nu- 
merical method. 



8. Migration 

It has become clear in recent years that protoplanets can 
be extremely mobile within their protoplanetary disk. Three 
types of migration can be distinguished: Type I, which con- 
cerns low-mass planets that do not open gaps in the disk , 
Type II for migration inside a clear gap, and Type III, for 
which the radial movement of the planet drives its migration 
dMass et & Panialoizou 2003). Because our planet stays at a 
fixed location, we are only dealing with Type I and Type II 
migration. 

For the torque calculations, we exclude material orbiting 
within the Roche lobe of the planet, and we assume an initial 
disk mass of 0.0035 Mq. 

8.1. Numerical Parameters 

Migration rates turn out to be very robust once the Roche 
lobe as well as the disk scale height is well resolved. Across 
the whole mass spectrum we found no significant torque dif- 
ferences for the various methods of source term integration 
and flux limiters. Even approximate extrapolation, which was 
shown to lead to spurious numerical evolution of the gap, does 
not affect the torque on the planet. This is because for low- 
mass planets approximate extrapolation does not alter the den- 
sity structure around the planet (see Sec. \1.2l . and for high- 
mass planets it only speeds up gap formation, which leads to 
Type II migration. 

As an example we show in Fig.[^the torque as a function 
of time for a 1/8 Mj planet for different flux limeters. All three 
limiters lead to a torque that is comparable to Type I migration. 
As with gap formation, the difference between p =1.5 andp = 
1.0 is larger than the difference between p — 1.5 and p = 2.0 
(see Fig.|5|l, but for the torque the effect is not significant. This 
shows that the sensitivity of the accretion rate on numerical 
parameters found in Sec. l7.3l is due to effects inside the Roche 
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Fig. 17. Migration time scales for p lanets of diff e rent m ass. 
Dashed line: analytical solution by iTanaka et al. ( 2002') for 
Type I migration. Dash-dotted line: result bv lWard (i 19971) for 
Type II migration. 



Lobe, while we explicitly excluded this region for the torque 
calculations of Fig. 1151 

This is further illustrated in Fig.^J where we again show 
the total torque on the planet, but now we exclude only the inner 
half of the Roche lobe. The first thing to note is that the minmod 
limiter now gives a significantly different result compared to 
the other two. We found similar behaviour for the accretion 
rate in Sec. 17. 31 



Secondly, we observe that the torque is less negative for 
p = 1.5 and p = 1 . 0. Thi s is the torque reversal also observed 
by Id] Angelo et a D (l2002h : material within the Roche lobe ex- 
erts a positive torque on the planet, slowing down migration. 
However, it is not clear exactly where the disk ends and the en- 
velope of the planet starts. A gas giant is usually assumed to 
fill its Roche lobe during formation (Pollack et al,.1996) , so all 
material orbiting inside the Roche lobe is part of the envelope. 
Therefore the question rises whether a planet should change 
its orbit due to its own envelope. On the other hand, impor- 
tant physical processes are not included in the cuiTent hydro- 
dynamical simulations of planet-disk interaction, namely heat 
transport (due to the isothermal assumption) and self-gravity. 
Both can have dramatic influence on the direct surroundings 
of the planet, where the density is highest and the temperature 
may differ significantly from the ambient disk, and therefore 
it is not feasible yet to make a connecti on between the one- 
dimensional planet formation models of IPoUack et alJ ( 1 19961) 
and the hydrodynamic models of planet-disk interaction. For 
now, we are only interested in comparing our values to previ- 
ous analytical and numerical results, which will be done in the 
next section. 



For the low-mass case of Type I migration there exists a nice 
analytical p rediction of the migr ation rate as a function of plan- 
etary mass jTanaka et al.ll2002h . For gap-opening planets, the 
migration rate should not de pend on planetary mass, only on 
the viscosity in the disk (Wajdl ll997h . Numerical simulations 
show that in between these two extremes the fastest migration 
takes place, in two-dimensional simulations Id' Angelo et alJ 
2002) as well as in three-dim ensional simulations ( Bate e t alJ 
(120031) : ^' Angelo et alJ j Eool v 

In Fig.^lwe show the migration time scale r,n defined by: 



(54) 



where the radial velocity of the planet due to a torque T is given 
by: 



2r. 



7^ 



(55) 



Here T is the torque due to the disk on the planet. 

From Fig.^]we see that we can reproduce Type I as well as 
Type II migration. The transition region extends approximately 
from Mp = 0.1 Mj to M„ = 1 .0 Mj, consistent with the three- 
dimensional results of bate et al. (2003). A planet of mass 0.1 
Mj has the fastest migration time scale of approximately 10"* 
years. 

Comparing Figs. 1141 and [TTl we see that for the low-mass 
planets the growth time scale is much shorter than the migra- 
tion time scale, while for the high-mass planets it is the other 
way around. In the intermediate case, both time scales are ap- 
proximately equal. 



9. Summary and conclusion 

We have presented a new method for modeling disk-planet in- 
teraction. Key features of the RODEO method are: it is a con- 
servative method, it treats shocks and discontinuities correctly, 
and it uses stationary extrapolation to integrate the geometrical 
and gravitational source terms. 

We found that the RODEO method performs very well on 
the problem of disk-planet interaction, not in the least because 
we do not experience s erious computational difficulties as for 
example in lKlevl(ll998h . We have shed some new light on this 
matter in that this instability can be seen to result from the 
failure of most hydrodynamic schemes to recognize stationary 
states. 

We find that the proces of gap formation crucially depends 
on source term integration. Only exact and hybrid extrapolation 
produce gaps of correct width. A diffusive flux limiter leads to 
a more pronounced inner edge of the gap. 

The accretion rates turn out to be very robust in the low- 
mass regime (Mp < 0.1 Mj) as well as in the high-mass regime 
(Mp >= 1.0 Mj): they are independent of numerical resolu- 
tion, accretion parameters and source term integration. A dif- 
ferent equation of state affects the accretion rates significantly: 
even a very low F of 1.001 that is used frequently to mimic 
isothermal flow the accretion rate drops by a factor of 2. 
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For intermediate mass planets the accretion rates are far less 
certain. They are dependent on the flux limiter used as well as 
on the way source terms are integrated. Within our method we 
find differences of 50%, and the differences bet ween our results 
and the similar study of D'Ange lo et alJ(l2002h are of the same 
order. This shows that the error bars on these values are very 
large. 

For the highest mass planets we find significantly lower ac- 
cretion rates than in previous numerical studies, limiting the 
maximum planet mass that can be reached to about 4 Mj. Note 
that this is in fact the first numerical study that consistently goes 
beyond this limiting mass so that no extrapolation is required. 

Migration rates are far more robust than accretion rates, as 
long as we limit the disk region that exerts a torque on the 
planet to outside the planetary Roche lobe. Therefore the re- 
gion where most differences arise in the accretion rates is lo- 
cated deep within the Roche lobe. 

We can nicely reproduce the analytical results for Type 1 
and Type II migration, with a transition region extending from 
Mp = 0.1 Mj to A/p = 1.0 Mj. In this transition region the 
fastest migration rate corresponds to a mass of 0.1 Mj. 

The RODEO method can also be used fo r two-fluid calcula- 
tions (see Paardekooper & Mellemal to model dust-gas 
interaction in protoplanetary disks. Also, the method can eas- 
ily be extended to three dimensions, and is also suited to treat 
an energy equation in a multi-dimensional set-up. We will con- 
sider these extensions in forthcoming papers. 
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Appendix A: Explicit expressions 

Here we give explicit expressions of the eigenvalues, eigenvec- 
tors and projection coefficients for the isothermal Euler equa- 
tions in cylindrical coordinates. 



A vector A = (A^, A^, A0) can be projected on these eigen- 
values using the following projection coefficients, found by 
solving the system Ca = A, where C is the matrix with the 
eigenvectors: 
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A.2. Azimuthal direction 

The Jacobian matrix A is given by: 
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A = I -VrV^, Vr 

i-vl 2v^^ 



The eigenvalues of this matrix are: 

Ai = v^- Cs/r 
A2 = U0 + Cs/r 

The corresponding eigenvectors: 

ei = {l,Vr,v^ - Cs/r) 

62 = [l,Vr,V^ + Cs/r) 

63 = (0,1,0) 



(A.5) 



(A.6) 



(A.7) 



A vector A = (Ap, A^, A^) can be projected on these eigen- 
values using the following projection coefficients, found by 
solving the system Ca = A, where C is the matrix with the 
eigenvectors: 
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A.1. Radial direction 

The Jacobian matrix A is given by: 
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-VrV^ Vtf, Vr 

The eigenvalues of this matrix are: 

Ai = Vr - Cs 
X2 — Vr + Cs 
A3 = Vr 

The corresponding eigenvectors: 

ei = {l,Vr - Cs,V^) 

62 (l,Wr +Cs,V^) 

63 = (0,0,1) 



(A.2) 



(A.3) 
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